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We performed molecular dynamics simulations of granular beds driven by a model hydrodynamic 
shear flow to elucidate general grain-scale mechanisms that determine the onset and cessation of 
sediment transport. By varying the Shields number (the nondimensional shear stress at the top 
of the bed) and particle Reynolds number (the ratio of particle inertia to viscous damping), we 
explore how variations of the fluid flow rate, particle inertia, and fluid viscosity affect the onset 
and cessation of bed motion. For low to moderate particle Reynolds numbers, a critical boundary 
separates mobile and static states. Transition times between these states diverge as this boundary is 
approached both from above and below. At high particle Reynolds number, inertial effects become 
dominant, and particle motion can be sustained well below flow rates at which mobilization of a 
static bed occurs. We also find that the onset of bed motion (for both low and high particle Reynolds 
numbers) is described by Weibullian weakest-link statistics, and thus is crucially dependent on the 
packing structure of the granular bed, even deep beneath the surface. 

PACS numbers: 92.40.Gc, 92.10.Wa, 45.70.Ht, 47.57.Gc 


I. INTRODUCTION 

Fluid flowing laterally over a granular bed exerts shear 
stress on the grains. This occurs in many natural settings 
and industrial ^plications, such as sediment transport 
in riverbeds [I|, Q and slurries pipes d, 13 . The ratio of 
the shear stress exerted by the fluid on the top of the 
bed to the buoyancy-corrected particle weight is known 
as the Shields number 0 Q. For small 0, no grain mo¬ 
tion occurs; at sufficiently large 0, however, grains can 
be entrained by the flow [g-fl^. Despite decades of re¬ 
search, the nature of the transition between static and 
mobilized granular beds is not well understood. The ge¬ 
ometric structure of the contact network in the bed deter¬ 
mines its mechanical strength (T6l - [l8l| - Bed mobilization 
is also strongly affected by the complex and unsteady 
fluid flow above the bed, as well as how strongly the fluid 
flow couples to the grains, as quantified by the particle 
Reynolds number Rep sEEini that measures how 
quickly grains equilibrate to the fluid flow. Weak stresses 
applied to the interior of the bed by fluid flowing through 
the pore spaces between grains may also play a role in bed 
mobilization (Til-I^ . Although empirical hydraulic mod¬ 
els capture some important aspects of sediment transport 
problems there is at present no fundamental under¬ 
standing of the relative contributions of these effects on 
the onset and cessation of grain motion. 

In this paper, we study a simplified model of a fluid- 
driven granular bed to clarify the essential physics at the 
onset of bed motion. In particular, we seek to understand 
the nature of the mobile-to-static and static-to-mobile 


transitions as a function of 0 and Rcp and to predict the 
parameter regime where hysteresis, defined as a finite dif¬ 
ference between 0o, above which a static bed will begin 
to move, and 0c, below which a mobile system will come 
to rest, occurs. 


We performed molecular dynamics (MD) simulations 
of a two-dimensional (2D) system composed of friction¬ 
less disks subjected to a simplified fluid flow that decays 
from a large value above the bed to a small value in¬ 
side the bed. Although our model is highly simplihed, 
with, for example, no explicit unsteadiness in the flow or 
friction between the grains, we find that 0c(i?ep) from 
the simulations is consistent with the behavior obtained 
from a l arg e collection of experiments on sediment trans¬ 
port @,[13, [HI- particular, we find plateau values 0c 
and 0c at low and high Rcp, with 0c > Q'f, and an 
intermediate Rcp regime that connects the two limiting 
values. In the low Rcp limit, there is a sharp transition 
at 0c between mobile and static beds in the infinite-time 
and infinite-system-size limits with no hysteresis. In the 
large Rep limit, we find significant hysteresis, since parti¬ 
cle inertia can sustain motion well below the 0o at which 
bed motion is initiated. We also find that the onset of bed 
motion at 0 > 0c for low Rep and 0 > 0o for high Rep 
depends strongly on system size and exhibits weakest- 
link statistics [Ij, Thus, the onset of bed motion in 
our system depends on the bed packing structure, even 
deep beneath the surface. 
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II. DETAILS OF THE MODEL 

We study a domain of width W that contains N/2 
large and N/2 small disks with diameter ratio 1.4. There 
is no upper boundary, and the lower boundary is rigid 
with infinite friction so that the horizontal velocities of 
all particles touching it are set to zero. We use periodic 
boundary conditions in the horizontal direction. The to¬ 
tal force on each particle is given by the vector sum of 
contact forces from other particles, a gravitational force, 
and a Stokes-drag-like force from a fluid that moves hor¬ 
izontally: 


niidi = '^F‘j - rmg'y + Bi[vo}{f)x - Vi], (1) 

i 

Here, rrii oc D/ is the particle mass, Di is the diameter of 
particle i, Vi and Si are the velocity and acceleration of 
particle z, rriig' is the buoyancy-corrected particle weight, 
Bi oc Di sets the drag on disk z, uq is a characteristic fluid 
velocity at the surface of a static bed, and f{r) is the 

fluid velocity at r. F// = K (l — ff fij is 

the pairwise repulsive contact force on disk i from disk j , 
where K is the particle stiffness, is the separation be¬ 
tween the centers of the particles, Dij = {Di -|- _Dj)/2, 
fij is the unit vector connecting their centers, and 9 
is the Heaviside step function. f{r), the fluid veloc¬ 
ity profile, varies smoothly from a large value above 
the bed to a small value inside the bed. We choose a 
form that depends only on the local packing fraction (j)i\ 
f{4>i) = e-K4'i-o.5] h controls the ratio of the mag¬ 

nitude of the fluid flow above and inside the bed. (pi is cal¬ 
culated in a small circular region with diameter Di -|- 2Di 
around each particle, as shown in Fig.[lja). We note that 
/ = 1 for = 0.5, a typical value at the bed surface. 
See Section IHI Al for force proHles in static and mobile 
states. 

Three nondimensional numbers govern the behavior of 
Eq. ([1]). We set the nondimensional stiffness > 3 x 
10^ to be sufficiently large that increasing it has no effect 
on our results. The other two nondimensional parameters 
can be written as 


0 = 


r = 


Bvq 

mg' 

B/m 

WTd' 


( 2 ) 

( 3 ) 


The Shields number 0 gives the dimensionless shear force 
at the top of a static bed. F is the ratio of the gravita¬ 
tional settling time Tg = \JD jg' to the viscous time scale 
m/B. Since the particle Reynolds number Rcp = 
where u is the kinematic viscosity, and the Stokes drag is 
proportional to pfi^D, where p/ is the fluid density, the 
ratio ^ oc ^Rcp {pg is the mass density of the 

grains) compares the inertia of grains entrained in the 
flow to the strength of the viscous drag. 


(a) 



(b) 



FIG. 1. (Color online) Layer-averaged fluid velocity Vf versus 
depth for static (a) and mobile (b) beds. Grains are subjected 
to a fluid drag force and a gravitational force —mg'y. The 
gray circle (left panel) defines the area used to calculate the 
local packing fraction (pi near the zth particle, which deter¬ 
mines the local fluid velocity. We also show the layer-averaged 
grain velocity Vg for a mobile bed in the right panel. 


To characterize flow onset and cessation in our system, 
we employed two protocols. In protocol A, to study the 
mobile-to-static transition, we distributed particles ran¬ 
domly throughout the domain and set a constant value 
of 0 for a total time of roughly lO^Tg. We consider the 
bed to be at rest when the maximum net particle acceler¬ 
ation Umax is below a threshold Othresh roughly one order 
of magnitude smaller than g' and roughly three orders of 
magnitude smaller than typical values for a moving bed. 
In protocol B, to understand the dynamics of the static- 
to-mobile transition, we begin with a static bed from pro¬ 
tocol A and slowly increase 0 in increments A0 = 0.010. 
If Omax < Othresh after roughly one inter-grain collision 
time, then 0 is increased. If Omax > Othresh, we keep 0 
constant until Omax < othresh- We designate a system as 
mobile under protocol B with slightly different criteria: 
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Omax/athresh > 10 End Vg > 0.0414, where Vg is the av¬ 
erage horizontal velocity of all grains and 14 = y/d'D is 
the settling velocity. These thresholds filter out small re¬ 
arrangement events, keeping only states with substantial 
grain motion. 


III. RESULTS AND DISCUSSION 

Figure [2] shows the boundaries between mobile and 
static beds as a function of 0 and Rep from simulations 
with 6 = 2, 4, and 6. We find a curve &c{Rep) above 
which the particles are unable to find a stable packing 
under protocol A. In the low-i?ep limit near 0c, grain 
motion is highly overdamped, and particles do not leave 
the bed. As Rcp is increased, the inertial effects of mo¬ 
bilized particles striking the bed make finding a stable 
configuration more difficult, decreasing 0c significantly. 
Figure [2] shows parameter values where grain motion does 
(blue circles) and does not (green squares) stop under 
protocol A. We then apply protocol B to stopped sys¬ 
tems and find another boundary Qo(Rep) that specifies 
when grain flow can be initiated. For low RCp, 0o < 0c, 
and bed motion initiated near 0o is temporary. Thus, 
at low Rcp, permanent grain motion is initiated at 0c in 
the large system limit, as we discuss below. However, in 
the high-i?ep limit where particle inertia is dominant, we 
observe significant hysteresis: grain motion is initiated at 
00 , which is well above the value of 0c where mobile par¬ 
ticles colliding with the bed can sustain bed motion. We 
also note that this result is consistent with where, 
in simulations of Aeolian transport at high i?ep, a sig¬ 
nificant perturbation or lift force was required near 0c 
to initiate grain motion. Temporary (filled black circles) 
and permanent (red crosses) motion under protocol B 
are also marked in Fig. [2] The basic nature of the flow 
diagram is insensitive to variations in 6, although the 
numerical values for 0o(i?ep) and 0c(i?ep) change. We 
note the similarities of 0c(i?ep) in Fig. [5] to the exper- 
imental and observational data compiled in 0,113, El, 
even though our model is highly simplified. Both display 
a plateau in the onset value of 0 at low Rcp, a decrease 
in the onset value at moderate Rcp, and a lower plateau 
value at high Rep. 

The b parameter sets the ratios between the fluid ve¬ 
locity above the bed (va), at the top of a static bed (vq), 
and in the bulk of the bed (vb) where grains are packed 
densely (at <j)i ~ 0.84). If a grain is well above the bed, 
(j)i Ri 0.1 (since the grain itself contributes to the local 
packing fraction), so the ratio of the fluid velocity for this 
grain to the fluid velocity at the top of an otherwise static 
bed (where (j)i ~ 0-5 and thus / = 1) is Va/v^ ~ 

Table U shows the ratios Va/vQ, vo/vb, and Va/vb for b = 2, 
4, and 6. Note that these ratios increase with b, and that 
varying b changes all ratios simultaneously. Obviously, 
a fluid flow model with an additional parameter could 
decouple the Va/vo and vo/vb ratios, but the model used 
here was chosen as a very simple way to interpolate be- 


(a) 





0/r2 Rep 


FIG. 2. (Color online) The flow diagrams. Shields number 
0 versus particle Reynolds number Rcp, with (a) b = 2, (b) 
& = 4, (c) and b = 6. Diagonal lines of data points correspond 
to lines of constant F. The symbols show systems that came 
to rest (o) or never stopped (□) under protocol A, and were 
permanently ( x ) or temporarily (•) mobile as 0 was increased 
under protocol B. The dashed line shows 0c, above which the 
inertial effects from particles entrained in the flow lead to 
sustained motion. The solid line indicates 0o, below which 
the system will never be mobilized. The two large black open 
circles mark the parameter values we study in Fig. [5] 
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TABLE I. The flow velocity ratios at different heights and the 
particle Reynolds number when ©o = 0c- Va, Vq, and vt, are 
the fluid velocities above the bed, at the surface of the bed, 
and in the interior of the bed, respectively. 


6 

Va/Vo 

Vo/Vb 

Va/Vb 

Rcp when 0o = 0c 

2 

2.23 

1.97 

4.39 

Rep « 50 

4 

4.95 

3.89 

19.3 

Rep « 10 

6 

11.02 

7.69 

84.8 

Rep « 3 


tween a large fluid velocity above the bed and a small 
fluid velocity inside the bed. Figure [2] shows that the 
flow diagram is qualitatively the same for 6 = 2, 4, and 
6. ©0 is relatively insensitive to b, but ©^ shifts to smaller 
Rcp and the gap between the plateau values at large and 
small Rcp widens with increasing b. 


A. Protocol A: Mobile-to-static transition 

Figure [3] shows data from protocol A near ©c, which 
we find to be nearly independent of system size, as we 
discuss below. In Fig. EJa), we show the time evolution 
of Vg multiplied by the fill height ND/W and normalized 
by Vg. At short times (red curve), Vg increases linearly 
with © and connects continuously to Vg = 0. This is 
the expected relation for frictionless granular systems: 
no motion occurs below the yield stress, while above the 
yield stress the strain rate increases as a power law in 
the difference between the applied and yield stresses [2^ . 
However, at long times we observe a discontinuity in Vg 
at ©c, as has also been observed in experiments [l3 as 
well as simulations of Aeolian transport [2^ and sheared 
frictional granular media [^ . The discontinuity in Vg 
moves toward ©c as time increases. The size of the dis¬ 
continuity and the slope Vg/Q both scale roughly linearly 
with Rcp, as shown in Fig-Efb). For Shields numbers be¬ 
low ©c, the grains settle into a stable packing in a time 
tg that diverges as © —>• ©c, as shown in Fig. He). 

Figure [1] shows that mobile grains are confined to a 
relatively small layer at the top of the bed. So, for fill 
heights studied here {ND/W > 8), we find the total flow 
rate VgND/W as a function of © and F is insensitive to 
the system size (recall, Vg is the average velocity of all 
grains). One might think that the discontinuity in the 
the flow rate near ©c may disappear in the large-system 
limit: small systems of grains are able to find a stable 
configuration, but very large systems will always have a 
weak spot that does not allow the system to stop. How¬ 
ever, this is not the case. Figure shows that the 
behavior near © = ©c is insensitive to the system size. 
The stopping time as a function of © is shown for three 
different system sizes, and the three curves are virtually 
identical. The mean and fluctuations of the stopping 
time both diverge as a power law as © —>■ ©c- The data 
shown is for F = 0.24 and 6 = 4, where ©c > ©o, and the 
power law exponent is roughly 2.5. As shown in Fig. He), 


(a) (b) 
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0 



(c) 
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0 

(d) 
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FIG. 3. (Color online) Data from protocol A near 0c. (a) 
Time evolution of the average grain velocity Vg{t) during pro¬ 
tocol A. The solid red curve characterizes grain motion at the 
time for grains to settle under no fluid flow tso , the black curve 
with open circles marks the end of the simulation, and dashed 
curves represents intermediate times (where red to black rep¬ 
resents increasing time). At long times, we observe a discon¬ 
tinuity in Vg with magnitude S and slope m at 0c. (b) J and 
m scale roughly linearly with Rep. (c) The time ts required 
for the grains to come to rest as a function of 0 for varying 
Rcp (left to right, blue to red, represents large to small Rep)\ 
symbols show F = 0.01 (4<), 0.02 (v), 0.05 (o), 0.1 (<), 0.15 
(*), 0.2 (>), 0.25 (□), 0.5 (A), and 1 (o). The lines show fits 
of (ts — tso) oc (0c — 0)“, where the values of a are marked 
next to each plot, (d) tg is plotted as a function 0 for 6 = 4, 
F = 0.24, and three system sizes. The symbols correspond 
to {W/D,N)'- black circles (100,800), blue squares (50,400), 
and red triangles (50,800). Data points show the mean of ten 
simulations, and error bars give the standard deviation. The 
inset shows a logarithmic plot of the mean of ts — tsO- The 
thick dashed line corresponds to (ts — tso) oc (0c — 0)~^ ®. 


this exponent is roughly 3 in the low Rcp limit, and it 
decreases with increasing Rcp to less than 1 in the high 
Rcp limit. ©0 is also shown, and it corresponds reason¬ 
ably well to where the divergence begins, which suggests 
that ©c — ©0 is related to the power law exponent of the 
divergence. 
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(a) 


(b) 


B. Protocol B: Static-to-mobile transition 




Fcimg 


FIG. 4. (Color online) Profiles of the layer-averaged instanta¬ 
neous fluid drag force (red lines) and the local pressure (black 
dashed lines) for b = 4 and F = 0.25, with (a) 0 = 0.2, (b) 
0.35, and (c) 0.4. 


In Figure m we show typical force profiles during pro¬ 
tocol A for grains under three different conditions. All 
panels have 6 = 4 and F = 0.25, but with varying 0. 
Solid red curves show the average force exerted by the 
fluid on the grains at a particular height. Black dashed 
curves show the average pressure due to grain-grain con¬ 
tacts as a function of height. To calculate the pressure, 
we first calculate the force moment tensor at par¬ 
ticle i by summing over the particles j that contact par¬ 
ticle i to obtain , where Ri is the 

particle radius, a and /3 represent components, and 
represents the /3 component of the branch vector connect¬ 
ing the center of particle i with the point of contact with 
particle j. The mean of the eigenvalues of this tensor pro¬ 
vides a grain-scale estimate of the local pressure, which 
we then average over the horizontal direction to obtain 
the average contact force Fc- We plot Fc in FigH] as a 
function of height from the lower boundary as a black 
dashed line, which has roughly slope 1 in all plots. Thus, 
the contact forces below the surface are dominated by 
the weight of grains above a particular layer. FigurelUJa) 
shows 0 = 0.2, which is below both 0o and 0c, so grains 
are not moving. Figure lU^b) shows 0 = 0.35, so the flow 
is metastable and grains will eventually come to a stop. 
Figure SKc) shows 0 = 0.4, so grains continue to move 
indefinitely. 


For 0 ^ 0c, a sufficiently large perturbation will lead 
to sustained bed motion. At high Rep, this motion never 
begins for 0 < 0o. But as the red crosses in Fig.[2]show, 
stable motion is not always initiated at these minimum 
values, but depends on the particular arrangement of 
grains in the bed. We show that sustained grain motion is 
consistent with Weibullian weakest link statistics [22|, [2^ 
and that failure events are always initiated at 0c for low 
Rep and 0o for high Rep for sufficiently large systems. 

First, we note that although the first large-scale motion 
of the grains always occurs at the top layer of grains, fail¬ 
ure events do not always originate there. If we measure 
the depth of the particle whose acceleration first exceeds 
^thresh at ©/, we find that it can occur at any depth yf 
below the bed surface, with a probability distribution pf 
that is proportional to the local applied fluid force as 
shown in Fig. [S] Essentially, the ratio of the probability 
of failure at the surface to the probability of failure in 
the bed is roughly given by the ratio between the fluid 
force at the top of the bed to that inside the bed vo/vg, 
as given in Table HI Thus, while large-scale particle mo¬ 
tion always begins at the surface of the bed, this motion 
is often correlated to small particle rearrangements that 
occur deep in the bed. 

Figure EKa) and (d) show that the probability distri¬ 
bution P{Of) approaches 6(0/ — 0c) for low Rep and 
6(0/ — 0o) for high Rep in the large system limit, where 
we varied both N and W to change the system size. The 
insets show that the distributions P{&f) collapse when 
rescaled by their mean values. The system-size depen¬ 
dence suggests a “weakest link” picture: at a given value 
of excess stress above 0c or 0o (for low and high Rep, 
respectively), there is a better chance of finding a suf¬ 
ficiently weak local grain arrangement somewhere in a 
larger system. If we consider the bed to be a composite 
system of M uncorrelated subsystems that fails if any of 
the subsystems fail, the cumulative distribution Cm{Q) 
for failure is related to that of a single subsystem (7(0) 
by 


1-Cm(0) = [1-^(0)]"^. 
If we assume a Weibull distribution [2^ [2^ 


(7(0) = 1 — exp 


/0-0C 

V p 



(4) 

(5) 


then Cm{F>) will have the same form with um = Oi and 
Pm = FigureEKa) and (d) show that P{Qf) = 

dC/dOf does indeed obey a Weibull distribution with 
shape parameter a ~ 2.6. 

Thus, if Eq. dS]) fits the data, and if a is constant but 
(3 oc as M is varied, then Eq. (SD applies, and 

global failure is caused by the failure of a single member 
of a collection of uncorrelated subsystems. Figure Etb) 
and (e) confirm this, showing that the mean Shields num¬ 
ber at flow onset 0/ scales as (0/ — 0c) oc for 
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(a) 



Vf /D 


FIG. 5. (Color online) The probability pf that bed failure 
occurs at depth yf (where yf = 0 is the top of the bed) is 
proportional to the local fluid force for (a) 6 = 2 and F = 0.1, 
(b) fo = 4 and F = 0.25, and (c) 6 = 6 and F = 0.5. The 
bed failure depth is determined as the depth of the particle 
whose acceleration first exceeds Uthresh (or a weighted aver¬ 
age of depths, if multiple grains meet this condition simulta¬ 
neously). The agreement between Pf{yj) and the local fluid 
force is good for b = A and 6 cases, but there is a deviation 
for the b — 2 case, such that Pf is smaller than expected near 
the lower boundary and bed surface. One explanation for this 
is that failure events are less localized, causing the weighted 
average of a; > Othresh to be more likely in the middle of the 
system. 


small Rcp and (Of — 0o) oc for large Ecp, where 

Meff = WesHes /Is the effective system size. This 
scaling means that larger systems are more likely to fail 
near 0c or 0o for small and large Rcp, respectively. How¬ 
ever, systems that fail near these minimum values are 
very slow to become fully mobilized. To demonstrate 
this, we consider the mobilization time defined as 
the time between the initial force imbalance that leads 
to failure and the time when Vg reaches its asymptotic 
value. As shown in Fig. m c) and (f), this time scale also 


diverges as 0 —>■ 0c for low Rcp and 0 Oq for high 
Rcp, and is independent of system size. 

To calculate Meg, we determine Wgg and Hgg as fol¬ 
lows. Since the vertical symmetry is broken by the fluid 
forcing profile, we calculate Ffeff by integrating the prob¬ 
ability of failure over the depth of the system, which is 
equal to the fluid force profile. That is, Hgg is the in¬ 
tegral of the prohles shown in Fig. [S] over the depth of 
the bed, since the relevant system size for the Weibullian 
weakest link scaling has to do with the probability of lo¬ 
cal failure, and, as previously mentioned, particles near 
the surface cause failure with a likelihood that is greater 
than particles beneath the surface by a factor vo/vb- We 
find that this form for Heg collapses the data for systems 
that are sufficiently large in the horizontal (periodic) di¬ 
rection. We also check the form of Heg for a case where 
6 = 2 and F = 0.1, such that vo/vb ~ 1.97, and we find 
it to collapse the data well according to the Weibullian 
scaling. The form of Hgg confirms our observation that 
surface grain motion can be initiated from deep beneath 
the surface, as it is also proportional to the local applied 
shear force. 

The horizontal dimension is symmetric, but we find 
finite system-size effects when W/D < where ^ is a 
horizontal correlation length that varies with Rep. As 
shown in Fig. [51 we find Meg = 0.91(0/ — 0c)““ for 
large systems, where Heg is calculated using the method 
described previously. However, when the horizontal di¬ 
mension becomes small, we observe that (0/ — 0c) is 
larger than expected, which corresponds to an effective 
system width Weg that is smaller than the real system 
width W. 

If Meff = 0.91(0/ - 0c)-“ = HegWeg/D^, and we 
assume that WegjW = (^(W/D), then we can write 

C(W/Z1)= 0.91(0/-0c)-“-^. (6) 

HegW 

A plot of this quantity is shown in Fig. [71 and the fit line 
corresponds to 


Q(W/D) = 1 - exp 


W 


(7) 


As Fig. 0 shows, the form WegjW = 1 — exp — ^ 
captures the finite size effects. We find a similar result 
for the the high Rep case in Fig. |5Kd)-(f), where F = 0.1 
and 6 = 4, but with a smaller value of ^ ~ 3. The form of 
Weg suggests a horizontal correlation length of roughly 
which is larger for low Rep (^ ~ 17) than for high Rcp 

(e^3). 


IV. SUMMARY 

In summary, we performed numerical simulations of a 
granular bed subjected to a simple fluid shear flow to 
understand general features of the initiation and cessa¬ 
tion of grain motion. The critical Shields number for the 
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FIG. 6. (Color online) Onset of bed motion is governed by Weibullian weakest-link statistics, (a)-(c) correspond to low Rcp 
with r = 0.25 and (d)-(f) correspond to high Rcp with F = 0.1. (a) and (d) show the probability distribution of the Shields 
number at bed failure 0/ for many different system sizes {W/D from 3.125 to 200, N from 25 to 1600, and WNjD from 8 
to 80). The dashed vertical lines define (a) 0c and (d) 0o. The insets show that P(0/) collapses when rescaled by 0/ — 0c, 
where 0/ is the mean of each distribution. The solid line is a Weibull distribution with shape parameter a ~ 2.6. (b) and 
(e) indicate that 0/ — 0c for low Rep and 0/ — 0o for high Rep scale with the effective system size . (c) and (f) show 

that tm, the mobilization time after bed failure, diverges near 0c and 0o, respectively, independent of system size. The insets 
show a logarithmic plot of tm. — tm,o versus (c) 0 — 0c and (f) 0 — 0o. The dashed lines show tm — tm,o oc (0 — 0c)(c) 
and tm — tm,o oc (0 — 0c)~° ® (f), and the thin vertical dashed line indicates 0c. Symbols *, A, o, <, >, □) correspond to 
different values of N (25, 50, 100, 200, 400, 800, and 1600, respectively) with varying WjD. Each data point represents an 
average of 20 simulations. 



W/D 


FIG. 7. (Color online) This plot shows the finite system size 
effects in the horizontal periodic direction (with 6 = 4 and 
r = 0.25, which are the same values for the low Rep data 
shown in Fig. [6Ka)-(c)). The horizontal axis is the system 
width in particle diameters, W/D. As we show in in Fig.[6Ka)- 
(c), we find Meg = 0.91(0/ — 0c)~“ for large systems. The 
vertical axis shows this quantity, 0.91(0/ — 0c)~“, divided by 
the product of the effective fill height Heg/D and the width 
in particle diameters W/D. This shows a good collapse, and 

the fit line corresponds to Wes/W = 1 — exp with 

^ = 16.7, and 95% confidence interval of roughly 15 < ^ < 20. 


onset of grain motion 0c(i?ep) from the simulations is 
consistent with the behavior from a large body of exper¬ 
imental results Q. At low i?ep, 0c(i?ep) separates mo¬ 
bile and static beds, but at high Rep, we observe signif¬ 
icant hysteresis as a consequence of particle inertia. We 
find that the onset of grain motion is directly connected 
to the packing structure, even deep in the bed where 
there is a weak but nonzero fluid stress [TqI - I^ . Our re¬ 
sults from this simple model clarify the essential physics 
governing the transition between mobile and static beds. 
In future work, additional effects such as turbulent flow, 
inter-grain friction, and nontrivial particle shape can be 
added one-by-one to determine their distinct effects on 
the onset and cessation of grain motion. 
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